{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise Set 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple Choice / True-False\n",
    "\n",
    "Correct answers are `marked like this`.\n",
    "\n",
    "**Q1.** When dealing with numerical instability, which of the following is NOT a recommended approach?\n",
    "\n",
    "> a. Using well-scaled and precise data\n",
    "\n",
    "> b. Ensuring matrix coefficients are within 6 orders of magnitude\n",
    "\n",
    "> c. `Increasing the Big M value as much as possible`\n",
    "\n",
    "> d. Setting the NumericFocus parameter to prioritize stability\n",
    "\n",
    "**Q2.** Which of the following statements about MIP Gap is accurate?\n",
    "\n",
    "> a. A lower MIP Gap always leads to significant solution improvements\n",
    "\n",
    "> b. `A higher MIP Gap can still yield a solution close to optimal, depending on the problem`\n",
    "\n",
    "> c. Setting MIP Gap to zero is ideal for all real-world applications\n",
    "\n",
    "> d. MIP Gap only affects runtime but not solution quality\n",
    "\n",
    "**Q3.** Which of the following is a potential benefit of using the Big M method in an optimization model?\n",
    "\n",
    "> a. `It helps establish conditional constraints based on binary variables`\n",
    "\n",
    "> b. It eliminates the need for slack variables\n",
    "\n",
    "> c. It reduces the numerical range required for coefficients\n",
    "\n",
    "> d. It ensures the model is solved to exact optimality\n",
    "\n",
    "**Q4.** An optimization model can sometimes be infeasible in exact arithmetic but feasible within the solver’s tolerances.\n",
    "\n",
    "> `True` \n",
    "\n",
    "> False\n",
    "\n",
    "**Q5.** For which of the following application areas is non-linear optimization used in practice?\n",
    "\n",
    "> a. Finance\n",
    "\n",
    "> b. Utilities\n",
    "\n",
    "> c. Engineering\n",
    "\n",
    "> d. `All of the above`\n",
    "\n",
    "**Q6.** Which type of optimization problem could have local optimal solutions that are not globally optimal?\n",
    "\n",
    "> a. Linear programs\n",
    "\n",
    "> b. Convex Quadratic Programs\n",
    "\n",
    "> c. `Non-convex Quadratic Programs`\n",
    "\n",
    "> d. None of the above\n",
    "\n",
    "**Q7.** Which of the following constraints defines a convex feasible region?\n",
    "\n",
    "> a. `Correct:` $y \\geq x^2$\n",
    "\n",
    "> b. $y \\leq x^2$\n",
    "\n",
    "> c. $y = x^2$\n",
    "\n",
    "> d. $y \\neq x^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple modeling / coding\n",
    "\n",
    "**Q8.** Show these constraints are either compatible (find a feasible solution) or incompatible (prove a contradiction) algebraically:\n",
    "> **Q8a.** \n",
    "\n",
    "\\begin{align*}\n",
    "x + y + z &\\ge 3   \\quad &(1)    \\\\\n",
    "2x + 3y - z &\\le 4   \\quad &(2)  \\\\\n",
    "x, y, z &\\ge 0\n",
    "\\end{align*}\n",
    "\n",
    "`Solution:`\n",
    "An example solution is to select $x = y = z = 1$.\n",
    "\n",
    "> **Q8b.** \n",
    "\n",
    "\\begin{align*}\n",
    "x_1 + x_2 + x_3 &= 6 \\quad &(1)      \\\\\n",
    "3x_1 + 2x_2 + x_3 &\\ge 20  \\quad &(2)\\\\\n",
    "x_1, x_2, x_3, &\\ge 0\n",
    "\\end{align*}\n",
    "\n",
    "`Solution:`\n",
    "\n",
    "One way to derive a contradiction:\n",
    "Multiply $(1)$ by 3, which gives $3x_1 + 3x_2 + 3x_3 = 18$. \n",
    "\n",
    "Since $x_1, x_2, x_3$ are all nonnegative,\n",
    "\\begin{align*}\n",
    "3x_1 + 3x_2 + 3x_3 \\ge 3x_1 + 2x_2 + x_3 \\ge 20\n",
    "\\end{align*}\n",
    "\n",
    "This reduces to:\n",
    "\n",
    "\\begin{align*}\n",
    "18 \\ge 3x_1 + &2x_2 + x_3 \\ge 20 \\\\\n",
    "18 &\\ge 20\n",
    "\\end{align*}\n",
    "\n",
    "Which is a contradiction. \n",
    "\n",
    "**Q9.** Translate the following into constraints using binary variables for deciding to invest in three stocks. Then determine if they are compatible as above. \n",
    "\n",
    "- I will invest in exactly two of the three stocks.\n",
    "- If I invest in stock A, then I will invest in stock B.\n",
    "- If I invest in stock B, then I will invest in stock C.\n",
    "\n",
    "`Solution:`\n",
    "Let $y_A, y_B, y_C \\in \\{0,1\\}$ be the three stocks.\n",
    "\\begin{align*}\n",
    "\\sum_{i \\in \\{A,B,C\\}} y_i &= 2 \\\\\n",
    "y_A &\\le y_B \\\\\n",
    "y_B &\\le y_C \\\\\n",
    "\\end{align*}\n",
    "$y_A = 0, y_B = 1, y_C = 1$ satisfies the constraints. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Back to the widget production problem\n",
    "Use the code below as the base model for this problem. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: gurobipy>=12.0.0 in /Users/yurchisin/opt/anaconda3/envs/gurobi_ml/lib/python3.11/site-packages (12.0.0)\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install 'gurobipy>=12.0.0'\n",
    "\n",
    "# Import packages\n",
    "import pandas as pd\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in transportation cost data\n",
    "path = 'https://raw.githubusercontent.com/Gurobi/modeling-examples/master/optimization101/Modeling_Session_1/'\n",
    "transp_cost = pd.read_csv(path + 'cost.csv')\n",
    "\n",
    "# get production and distribution locations from data frame\n",
    "production = list(transp_cost['production'].unique())\n",
    "distribution = list(transp_cost['distribution'].unique())\n",
    "transp_cost = transp_cost.set_index(['production','distribution']).squeeze()\n",
    "\n",
    "max_prod = pd.Series([180,200,140,80,180], index = production, name = \"max_production\")\n",
    "n_demand = pd.Series([89,95,121,101,116,181], index = distribution, name = \"demand\")\n",
    "\n",
    "# the min production is a fraction of the max\n",
    "frac = 0.75\n",
    "C = 30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Restricted license - for non-production use only - expires 2026-11-23\n",
      "The original model had a total cost of 1748.42\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>cost</th>\n",
       "      <th>shipment</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>production</th>\n",
       "      <th>distribution</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Baltimore</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>5.96</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Richmond</th>\n",
       "      <td>1.96</td>\n",
       "      <td>116.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Cleveland</th>\n",
       "      <th>Columbia</th>\n",
       "      <td>2.43</td>\n",
       "      <td>89.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>2.37</td>\n",
       "      <td>95.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <th>St. Louis</th>\n",
       "      <td>2.92</td>\n",
       "      <td>140.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <th>Nashville</th>\n",
       "      <td>1.53</td>\n",
       "      <td>71.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Charleston</th>\n",
       "      <th>Lexington</th>\n",
       "      <td>1.61</td>\n",
       "      <td>121.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>St. Louis</th>\n",
       "      <td>4.60</td>\n",
       "      <td>41.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          cost  shipment\n",
       "production  distribution                \n",
       "Baltimore   Nashville     5.96      30.0\n",
       "            Richmond      1.96     116.0\n",
       "Cleveland   Columbia      2.43      89.0\n",
       "            Indianapolis  2.37      95.0\n",
       "Little Rock St. Louis     2.92     140.0\n",
       "Birmingham  Nashville     1.53      71.0\n",
       "Charleston  Lexington     1.61     121.0\n",
       "            St. Louis     4.60      41.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# gurobipy code for model\n",
    "m = gp.Model('widgets')\n",
    "m.setParam('OutputFlag',0)\n",
    "\n",
    "# decision variables - reminder that we ended on x being a semi-continuous variable type\n",
    "x = m.addVars(production, distribution, vtype=GRB.SEMICONT, lb = C, name = 'prod_ship')\n",
    "\n",
    "# constraints\n",
    "can_produce = m.addConstrs((gp.quicksum(x[p,d] for d in distribution) <= max_prod[p] for p in production), name = 'can_produce')\n",
    "must_produce = m.addConstrs((gp.quicksum(x[p,d] for d in distribution) >= frac*max_prod[p] for p in production), name = 'must_produce')\n",
    "meet_demand = m.addConstrs(x.sum('*', d) >= n_demand[d] for d in distribution)\n",
    "\n",
    "#objective\n",
    "m.setObjective(gp.quicksum(transp_cost[p,d]*x[p,d] for p in production for d in distribution), GRB.MINIMIZE)\n",
    "\n",
    "m.optimize()\n",
    "\n",
    "x_values = pd.Series(m.getAttr('X', x), name = \"shipment\", index = transp_cost.index)\n",
    "sol = pd.concat([transp_cost, x_values], axis=1)\n",
    "\n",
    "print(f\"The original model had a total cost of {round(m.ObjVal,2)}\")\n",
    "sol[sol.shipment > 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear optimization modeling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quadratic expression\n",
    "\n",
    "In order to produce more widgets, production facilities must run their machines at faster speeds. This consumes more energy and produces more waste at a level that grows with each widget produced. Running their machines costs a production facility $0.01 times the square of the number of widgets produced. Modify the objective function to include this additional production cost. Here is the some code to get you started - it includes the original transportation cost expression that you will be adding to:\n",
    "\n",
    "```python\n",
    "m.setObjective(gp.quicksum(transp_cost[p,d]*x[p,d] for p in production for d in distribution) + ..., GRB.MINIMIZE)\n",
    "```\n",
    "Solve the model with the updated objective and see how the cost has changed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The updated model had a total cost of 2808.99\n"
     ]
    }
   ],
   "source": [
    "# One solution\n",
    "m.setObjective(gp.quicksum(transp_cost[p,d]*x[p,d] for p in production for d in distribution) + 0.01 * gp.quicksum(x.sum(p, '*') * x.sum(p, '*') for p in production), GRB.MINIMIZE)\n",
    "m.optimize()\n",
    "print(f\"The updated model had a total cost of {round(m.ObjVal,2)}\")\n",
    "\n",
    "# Another solution\n",
    "# Introduces auxiliary variables representing total production for each facility\n",
    "# produced = m.addVars(production, name='total_produced')\n",
    "# m.addConstrs((produced[p] == x.sum(p, '*') for p in production), name='set_produced')\n",
    "# m.setObjective(gp.quicksum(transp_cost[p,d]*x[p,d] for p in production for d in distribution) + 0.01 * gp.quicksum(produced[p] * produced[p] for p in production), GRB.MINIMIZE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nonlinear expression 1\n",
    "\n",
    "The widgets will now be shipped using hyper-efficient vehicles. The cost to ship $x_{ij}$ widgets from production facility $i$ to distribution center $j$ is now equal to the corresponding transportation cost times $\\log(1 + x_{ij})$. Modify the objective function to be equal to this new transportation cost. Some code to get you started:\n",
    "\n",
    "```python\n",
    "prod_cost = m.addVar(name='prod_cost')\n",
    "m.addConstr(prod_cost ==..., name='set_prod_cost')\n",
    "m.setObjective(prod_cost, GRB.MINIMIZE)\n",
    "```\n",
    "Again, resolve the model and see how the cost has changed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The updated model had a total cost of 97.47\n"
     ]
    }
   ],
   "source": [
    "# One solution\n",
    "prod_cost = m.addVar(name='prod_cost')\n",
    "m.addConstr(prod_cost == gp.quicksum(transp_cost[p,d] * gp.nlfunc.log(1 + x[p, d]) for p in production for d in distribution), name='set_prod_cost')\n",
    "m.setObjective(prod_cost, GRB.MINIMIZE)\n",
    "m.optimize()\n",
    "print(f\"The updated model had a total cost of {round(m.ObjVal,2)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nonlinear expression 2\n",
    "\n",
    "A very strange company rule stipulates that the following mathematical relationship must hold:\n",
    "\n",
    "$ \\log\\Large(\\frac{4 \\cdot x_{\\textrm{Baltimore,Indianapolis}}}{\\sum_d x_{\\textrm{Baltimore},d}})\\normalsize \\geq 0$\n",
    "\n",
    "Add this constraint to the model using a nonlinear expression.\n",
    "\n",
    "Hint: Start by adding a non-negative slack variable $s$ and writing the constraint as $s = \\log\\Large(\\frac{4 \\cdot x_{\\textrm{Baltimore,Indianapolis}}}{\\sum_d x_{\\textrm{Baltimore},d}})$\n",
    "\n",
    "Solve the model and see how the cost has changed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The updated model had a total cost of 106.92\n"
     ]
    }
   ],
   "source": [
    "# One solution\n",
    "s = m.addVar(name=\"s\")\n",
    "NL2 = m.addConstr(s == gp.nlfunc.log(4 * x['Baltimore', 'Indianapolis'] / x.sum('Baltimore', '*')), name='strange_company_rule')\n",
    "m.optimize()\n",
    "print(f\"The updated model had a total cost of {round(m.ObjVal,2)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### Linearizing a nonlinear expression\n",
    "\n",
    "Wait a minute...is the nonlinear constraint we just added unnecessarily complicated? Rewrite the nonlinear constraint using only linear expressions. Describe the purpose of the constraint in words.\n",
    "\n",
    "Solve the model again and inspect the solution. Verify by inspection that it satisfies this newly added constraint.\n",
    "\n",
    "Using the properties of logarithms, the constraint can be rewritten as:\n",
    "$ x_{\\textrm{Baltimore,Indianapolis}}\\geq \\frac 14\\sum_d x_{\\textrm{Baltimore},d}$\n",
    "\n",
    "\\begin{align*}\n",
    "\\log(\\frac{4 \\cdot x_{\\textrm{Baltimore,Indianapolis}}}{\\sum_d x_{\\textrm{Baltimore},d}})\\normalsize &\\geq 0 \\\\\n",
    "\\frac{4 \\cdot x_{\\textrm{Baltimore,Indianapolis}}}{\\sum_d x_{\\textrm{Baltimore},d}} &\\geq 1 \\\\\n",
    "x_{\\textrm{Baltimore,Indianapolis}}&\\geq \\frac 14\\sum_d x_{\\textrm{Baltimore},d}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The updated model had a total cost of 106.92\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>cost</th>\n",
       "      <th>shipment</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>production</th>\n",
       "      <th>distribution</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Baltimore</th>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>5.09</td>\n",
       "      <td>38.666667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Richmond</th>\n",
       "      <td>1.96</td>\n",
       "      <td>116.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Cleveland</th>\n",
       "      <th>Columbia</th>\n",
       "      <td>2.43</td>\n",
       "      <td>89.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Nashville</th>\n",
       "      <td>4.13</td>\n",
       "      <td>101.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Little Rock</th>\n",
       "      <th>St. Louis</th>\n",
       "      <td>2.92</td>\n",
       "      <td>121.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Birmingham</th>\n",
       "      <th>St. Louis</th>\n",
       "      <td>4.01</td>\n",
       "      <td>60.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">Charleston</th>\n",
       "      <th>Indianapolis</th>\n",
       "      <td>2.61</td>\n",
       "      <td>56.333333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lexington</th>\n",
       "      <td>1.61</td>\n",
       "      <td>121.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          cost    shipment\n",
       "production  distribution                  \n",
       "Baltimore   Indianapolis  5.09   38.666667\n",
       "            Richmond      1.96  116.000000\n",
       "Cleveland   Columbia      2.43   89.000000\n",
       "            Nashville     4.13  101.000000\n",
       "Little Rock St. Louis     2.92  121.000000\n",
       "Birmingham  St. Louis     4.01   60.000000\n",
       "Charleston  Indianapolis  2.61   56.333333\n",
       "            Lexington     1.61  121.000000"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# One solution\n",
    "# First, remove the complicated form of the constraint\n",
    "m.remove(NL2)\n",
    "\n",
    "# Meaning of the constraint: at least 25% of Baltimore's total production must be sent to Indianapolis\n",
    "m.addConstr(x['Baltimore', 'Indianapolis'] >= 0.25 * x.sum('Baltimore', '*'), name='strange_company_rule_linearized')\n",
    "m.optimize()\n",
    "x_values = pd.Series(m.getAttr('X', x), name = \"shipment\", index = transp_cost.index)\n",
    "sol = pd.concat([transp_cost, x_values], axis=1)\n",
    "\n",
    "print(f\"The updated model had a total cost of {round(m.ObjVal,2)}\")\n",
    "sol[sol.shipment > 0]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gurobi_ml",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
